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INTRODUCTION 


Traditionally,  Lagrangian  codes  have  been  used  to  simulate  material  response  when  the 
amount  of  deformation  is  small.  When  the  deformation  is  large,  Eulerian  calculations  have 
been  employed.  The  Lagrangian  calculation  is  more  accurate  -  the  Eulerian  calculation  has 
greater  applicability.  These  strengths  and  weaknesses  are  due  to  the  convective  derivative 
which  is  absent  in  the  equations  written  in  the  moving  Lagrangian  frame.  Numerical 
treatment  of  this  advection  term  is  difficult  and  introduces  inaccuracies  into  the  calculation. 
However,  if  these  errors  can  be  made  small,  the  Eulerian  calculation  can  be  used  to  treat  a 
variety  of  high  strain  phenomena. 

Various  methods  have  been  devised  in  order  to  achieve  the  best  features  of  both  approaches. 
These  “hybrid”  techniques  normally  use  two  grids,  one  Lagrangian  -  the  other  Eulerian, 
with  information  exchanged  between  them.  These  mappings  add  a  good  deal  of  complexity 
to  the  calculation  and  can  also  introduce  inaccuracies.  Nevertheless,  many  hybrid 
techniques  have  been  successful  and  are  widely  used  today 

Unique  in  computational  fluid  dynamics  is  Smoothed  Particle  Hydrodynamics  (SPH).  This 
technique  uses  no  underlying  grid  -  it  is  a  pure  Lagrangian  panicle  method  invented  by 
Lucy  [1],  Gingold  [2,3],  Monaghan  [4,5,6]  and  Benz  [7],  The  absence  of  a  mesh  means  that 
large  deformations  can  be  computed  in  a  pure  Lagrangian  frame.  It  is  for  this  reason  that 
SPH  has  the  potential  to  be  a  valuable  computational  tool.  Although  SPH  has  been  proven 
an  excellent  computational  tool  for  astrophysical  applications,  its  ability  to  treat  typical 
hydrocode  production  problems  is  largely  untested  at  this  point.  The  method  is  just  nosv 


being  applied  to  a  broad  range  of  problems  where  its  strengths  and  weaknesses  are  sure  to _ j  1 
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THE  METHOD 


The  foundation  of  Smoothed  Particle  Hydrodynamics  is  interpolation  theory.  The 
conservation  laws  of  continuum  fluid  dynamics,  in  the  form  of  partial  differential 
equations,  are  transformed  into  integral  equations  through  the  use  of  an  interpolation 
function  that  gives  the  “kernel  estimate”  of  the  field  variables  at  a  point.  Computationally, 
information  is  known  only  at  discrete  points,  so  that  the  integrals  are  evaluated  as  sums  over 
neighboring  points.  These  “interactions"  result  in  a  net  force  which  will  accelerate  the 
“particle”.  The  reason  that  an  underlying  grid  is  not  needed  is  that  functions  are  evaluated 
using  their  value  at  the  discrete  points  (particles)  and  an  interpolation  kernel.  An 
integration  by  parts  then  moves  spatial  derivatives  from  operating  on  the  physical  quantities 
to  operating  on  the  interpolation  kernel  which  is  an  analytic.  These  concepts  will  now  be 
described  more  fully.  Consider  a  function /,  a  kernel  W  which  has  a  width  measured  by  the 
parameter  h,  and  the  following  equation: 


<m  >=  j^r-r'.AMr'yr'  .  (1) 

If  the  integral  of  W  is  normalized  to  unity,  then  it  follows  that 

<  M  >  ~o  ~  Kr)  ■  (2) 

Relation  (1)  therefore  defines  the  kernel  estimate  <f>  of  /.  If  W  is  the  Dirac  delta  function 
then  we  have  the  equality  </>=/.  Now  suppose  that /is  known  only  at  N  discrete  points  that 
are  spatially  distributed  according  to  the  number  density  distribution: 

N 

«(«•)  =  Z<5(r_r')  (3) 

/-I 

If  we  associate  with  particle  j  a  volume 
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thus  introducing  the  concept  of  particle  mass  (m),  it  follows  from  (1)  that 


<M>  =  X  fi  Wfr-ry.A)—  (5) 

i  Qj 

This  equation  defines  a  procedure  for  transforming  integral  equations  to  particle  equations 
and  is  therefore  called  “integral  evaluation  by  the  particle  method.”  A  detailed  discussion 
of  the  theory  of  SPH  is  given  by  Benz  [7]. 


DERIVATION  OF  THE  SPH  EQUATIONS 

The  conservation  equations  of  continuum  mechanics  are: 


and 


dQ  dlf 

~dt  =  ~9  ~d? 


(6) 


dU° 

1 

do “P 

V) 

dt 

Q 

dxP 

dE_ 

dip 

(8) 

dt 

Q 

dx? 

dxa 

9* 

=  Ua 

(9) 

dt 


Dependent  variables  are  the  scalar  density  (p)  and  specific  internal  energy  (E),  the  velocity 

vector  velocity  Ua ,  and  the  stress  tensor  .  The  independent  variables  are  the  spatial 
coordinates  (*)  and  the  time  ( t ),  and  the  total  time  derivative  ( d/dt )  is  taken  in  the  moving 
Lagrangian  frame.  Summation  over  repeated  greek  indicies  is  implied.  Let  us  now  cast 
equations  (6-8)  into  the  SPH  framework  by  applying  the  procedure  outlined  above.  First, 
rewrite  the  momentum  and  energy  equations,  so  that  the  density  (p)  appears  inside  the  spatial 
derivative  operator,  then  find  the  kernel  estimate.  The  result  is 
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We  now  linearize  these  equations  by  taking  integrals  of  products  equal  to  products  of 
integrals  (a  second  order  accurate  approximation),  giving 
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The  right-hand-sides  of  these  equations  are  now  integrated  by  parts,  assuming  W 
approaches  zero  fast  enough  that  the  surface  terms  vanish. 
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Finally,  the  integrals  are  evaluated  by  the  particle  method,  Eq.  (5),  to  give 


dQi 

dt 

(19) 

dUf 
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(20) 

and 

dEi 

dt 

=  ■ 
t'r  , 

(21) 

We  have  introduced  the  notation  and  dWjj/dx f  =  In  obtaining  (19) 

we  subtracted  from  (13)  the  following  term, 


which  is  zero  because  the  kernel  vanishes  at  intinity.  In  this  way  we  introduce  velocity 
differences  into  the  density  calculation,  which  is  desirable  and  consistant  with  the  energy 
calculation  in  (21).  Eqs.  (19-21)  are  the  conservation  laws  of  continuum  dynamics  written 
in  the  SPH  framework.  A  given  particle  i  has  a  density  determined  by  (19),  an  acceleration 
obtained  from  (20)  and  an  internal  energy  change  given  by  (21).  The  summations  are  over 
neighboring  j  particles.  These  equations  are  not  unique.  Several  other  forms  of  particle 
equations  can  be  derived  using  various  mathematical  manipulations.  Some  of  these  are 
discussed  by  Monaghan  (8). 


THE  DENSITY  CALCULATION 

It  is  important  to  recognize  that  (19)  is  not  the  density  calculation  that  normally  appears  in 
the  SPH  literature.  It  is  more  in  the  spirit  of  Smoothed  Particle  Hydrodynamics  to  compute 
the  density  using  the  equation  obtained  by  substituting  p  for  </>  in  (5).  namely 
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(23) 


Qi  = 

j 

With  this  equation  only  particle  coordinates  and  masses  are  required  to  compute  the  density, 
and  the  continuity  equation  (6)  is  automatically  satisfied.  The  disadvantage  of  using  (23)  is 
edge  effects  -  particles  near  a  free  surface  appear  underdense  and  therefore  in  tension, 
causing  motion.  Benz  [7]  discusses  several  possible  solutions  to  this  problem  including 
spacing  modification,  ghost  particles,  and  initial  relaxation  and  the  use  of  (19).  It  is  worthy 
to  note  that  differentiation  of  (23)  leads  to 


which  differs  from  (19)  only  in  that  Qi  appears  in  the  denominator  rather  than  Qj.  We  have 
not  yet  explored  the  consequences  of  using  (24)  in  place  of  ( 1 9).  The  difference  is  of  the  same 
order  as  the  difference  between  the  product  of  the  expected  values  and  the  expected  value  of 
the  product. 

ARTIFICIAL  VISCOSITY  &  WALL  HEATING 

As  they  stand,  Eqs.  (19-21)  yield  large  unphysical  oscillations  near  shocks.  In  fact,  any 
numerical  solution  of  the  continuum  equations  will  exibit  this  behavior  because  tne 
dissipative  terms  have  been  omitted.  Variations  of  physical  quantities  across  shocks  in 
nature  are  far  to  sharp  to  be  captured  by  numerical  techniques.  Von  Neumann  and 
Richtmyer  [9]  invented  “artifical  viscosity"  which  acts  to  smooth  shocks  over  a  few 
resolution  lengths  and  stabilize  numerical  solutions.  The  additional  term  is  introduced  into 
the  equations  as  an  artifical  viscous  pressure  IT  We  follow  Monaghan  and  Gingold  [4]  who 
derived  the  following  artificial  viscous  pressure  for  SPH: 
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if  (u,  -  UjYir,  -  Fj)  <  0 


(25) 


where 


and 


n,  -  I 
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QH 
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otherwise 


h(ui-u,)'(fi-rj ) 


Qj  =  (C/  4-  Cy)/2  £>i)  =  (g»,  +  g»y)/2 


(26) 


(27) 


The  three  parameters  appearing  in  these  equations  have  typtica!  values  of  o:=0.5  ,  p=1.0  and 
e=0.1.  The  linear  term  in  (25)  uses  the  sound  speed  (c).  This  artifical  viscosity  gives 
satisfactory  results  in  most  cases,  but  under  some  severe  conditions  it  fails  to  remove 
spurious  heating.  An  example  of  this  is  when  a  stream  of  gas  is  brought  to  rest  against  a 
rigid  wall.  Noh  [10]  was  able  to  improve  numerical  solutions  in  such  cases  dramatically  by 
adding  an  artificial  heat  conduction  term  to  the  energy  equation.  Monaghan  [11]  derived 
the  SPH  analog  of  Noh’s  “wall  heating”  term  in  which  the  net  artifical  heat  flux  at  particle  i 
is  given  by 


(28) 

where 

£  =  gihc  +  g^lV-vl-V-v)  ,  (29) 

and 

&-(&  +  &)/: 2  ea  =  (<?,-  +  a)/2  0>o) 

Suitable  values  of  the  two  parameters  appearing  in  (29)  are  g i  =  0.5  and  gi  —  10. 
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CONSTITUTIVE  RELATIONS 


The  stress  tensor  apperaring  in  Eqs.  (20)  and  (21)  is  defined  in  terms  of  an  isotropic  part 
which  is  the  pressure  (P)  and  the  traceless  symmetric  deviatoric  stress  (S): 

=  Pft#  -  S .  (31) 


The  pressure  is  normally  computed  using  an  equation  of  state  having  functional  form 
P=P(p,E),  such  as  the  Mie-Gruneisen  equation  for  solids  and  gamma-law  for  gases. 


Mie-Gruneisen 


P  =  P^l-r^  +  T’QE 


Ph 


aQt 7  +  b0r?  +  cQr? 
a0V 


t]  >  0 
t]  <  0 


(32) 


Ideal  Gas  P  =  (y-  l)p£  (33) 

The  subscript  ”//”  refers  to  the  Hugoniot  curve,  while  Ti=p/po~l  is  used  to  represent  the 
compression  and  r’=Tp/po.  For  the  anisotropic  part  of  (31)  we  write  a  prognostic  equation  for 
the  deviatoric  stress  assuming  small  displacements 

Safi  =  ^  (34) 

where  p.  is  the  shear  modulus  and  £  is  the  traceless  rate  of  strain.  However,  for  finite 
displacements  this  equation  is  not  material  frame  indifferent  [12]  ,  that  is,  the  material 
response  will  depend  in  an  unphysicai  way  on  rotations  (and  possibly  translations)  of  the 
material  and  of  the  observer  describing  it.  A  variety  of  frame  indifferent  stress  rates  have 
been  formulated.  Herrmann  [13]  examines  the  relative  merits  of  several  of  these.  The 
Jaumann  rate  is  the  most  widely  used  in  codes  and  we  adopt  it  also.  With  the  Jaumann  rate, 
our  constitutive  equation  is 
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Safi  _  yyRPy  _  ^ay  _  ^  (35) 

The  strain  rate  and  rotation  rate  tensors  that  have  been  used  are  defined  as  follows: 


2  \  dx&  dxc 


/ 


2  \  fa?  dxa 


(36) 


Particle  equations  for  (35)  are  obtained  by  Libersky  and  Petschek  [14]  in  a  manner  similar 
to  that  of  (19,20,21) 


^  -  ■sffif'  -  if  Ftp (Uf-tfWM+ty-tfww-jDiSI* 


(37) 


The  divergence  is  already  determined  by  (19),  A  =  ~Qi/Qi  ,  and  the  rotation  rate  is 


The  plastic  flow  regime  is  determined  by  the  von-Mieses  criterion  when  the  second  stress 


invariant  J 2  =  S^S0^  exceeds  the  known  flow  stess  (Yo).  The  individual  deviators  are  then 
brought  back  to  the  flow  surface. 


Sa$  ® 


(39) 


A  more  accurate  treatment  f^r  most  metals,  not  yet  implemented  in  our  code,  is  obtained  by 
computing  a  history  sensitive  flow  stress,  rather  than  a  predetermined  fixed  value  described 
above.  The  Johnson-Cook  model  [15],  for  example,  takes  into  account  thermal  softening, 
strain  harding  and  strain  rate  effects  on  the  equivalent  flow  stress.  This  more  sophisticated 
model  contains  seven  strength  related  parameters.  The  elastic-perfectly  plastic  constitutive 
model  described  above  contains  two  parameters,  the  shear  modulus  (jx)and  the  plastic  yield 
stress  (Yo). 
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THE  KERNEL 


The  interpolation  kernel  or  smoothing  function  most  widely  used  in  SPH  is  the  B-spline  WV 


/ 


W<(v,h)  = 


M 


14 


0  <  v  <  1 

1  <  v  <  2 


0 


otherwise 


v  =  \ri~rj\/h. 


(40) 


The  fractional  coefficients  appearing  in  (40)  assure  proper  normalization  and  continuity. 
This  kernel  interpolates  to  second  order  in  h  and  is  always  positive  in  the  range  of  interest. 
The  kernel  also  has  compact  support,  that  is,  it  goes  to  zero  at  a  distance  2h  from  its  peak. 
This  provides  a  clear  limit  on  rhe  number  of  neighbor  particles.  A  Gaussian  kernel  is  second 
order  accurate  and  positive  definite,  but  the  lack  of  compact  support  necessitates  an 
artificial  cut-off,  often  taken  at  v=3,  making  it  less  a  less  efficient  choice.  Higher  order 
interpolation  kernels  exist  [5]  but  are  not  always  positive  definite. 


TIME  INTEGRATION 

Eqs.  (19,20,21,37)  are  integrated  using  a  standard  leap-frog  algoritihm  (16]  with  time  step 
8 1,  calculated  from  the  configuration  at  time  t,  tc  advance  the  field  variabled  to  1 4-  8/.  We  will 
switch  from  superscript  tensor  indicies  to  subscripts  here  in  order  to  accomodate  the  standard 
superscript  representation  of  the  time  stepping  in  which  n  indicates  the  current  time  t  and 
n  +  1  indicates  the  advanced  time  /  +  8t. 
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0*"=e*(l-D6t) 

(41) 

(fafVl  =  ITa*  +  V2{6tn  +  6tn-')F 

(42) 

Entl  =  £"  +dtnG 

(43) 

stf-s^  +  drH 

(44) 

xXx  =xna  +  Una"*6tn 

(45) 

In  these  equations  F,G  and  H  represent  the  total  acceleration  of,  work  per  unit  mass  and 
stress  rate  on  a  particle  as  determined  by  the  interactions  with  neighbor  particles.  The 
accuracy  of  the  leap-frog  scheme  is  second  order  in  time  and  its  stability  is  guaranteed  by 
using  the  CFL  condition  to  determine  the  size  of  the  time  step  fir.  We  find  the  minimum  over 
all  particles  of  c oh/(c  +s),  where  c  is  the  adiabatic  sound  speed,  s  is  the  particle  speed,  h  is  the 
smoothing  length  and  u  a  constant  factor.  Choosing  u)  =  0.3  seems  adequate. 

CODE  ARCHITECTURE 

MAGI  differs  from  most  codes  in  that  it  was  designed  from  the  beginning  for  application  to 
very  large  problems  on  vecto»  supercomputers.  Strategies  for  the  efficient  implementation 
of  the  SPH  method  were  considered  and  implemented  throughout  the  design  and  coding  of 
MAGI.  Vectorization,  an  efficient  neighbor  searcher,  accomodation  for  the  symmetry  of  the 
particle  interactions,  and  activity  flags  were  all  exploited  for  efficiency  and  reduced 
computation. 

Activity  flags,  which  mark  particles  experiencing  motion  cr  acceleration  are  used  to  gam 
efficiency.  Only  active  particles  and  those  within  a  two  smoothing  lengths  of  active  particles 
need  to  be  updated.  Stationary,  unshocked  material  remains  inactive  until  impacted  by 
moving  material  or  accelerated  by  a  shock  or  stress  wave.  This  capability  results  in 
significant  savings  in  computer  time  for  problems  that  contain  a  large  number  of  panicles 
that  are  initially  unaffected  by  the  impact. 


11 


MAGI  consists  of  a  group  of  subroutines,  partitioned  by  task,  that  are  called  in  logical 
sequence  each  computational  cycle.  They  accomplish  the  following  tasks  that  are  basic  to 
the  SPH  algorithm.  (1)  Compute  the  panicle  interaction  sums  on  each  panicle  to  determine 
the  accelerations,  strain  rates,  and  energy  increment.  (2)  Update  the  velocities,  energies, 
stresses,  and  density  on  each  panicle.  (3)  Update  the  pressure  on  each  panicle  using  the 
new  density  and  energy.  (4)  Advance  the  particle  positions.  The  subroutines  that  perform 
most  tasks  consist  of  a  single  FORTRAN  DO  loop  that  is  easily  vectorized.  The  interaction 
subroutine,  however,  is  an  exception.  It  consists  of  nested  loops,  an  outer  unvectorized  loop 
over  all  active  panicles,  and  two  inner  loops  that  are  indexed  over  the  list  of  neighbors  that  is 
returned  by  the  linked  list  described  below.  Each  of  the  inner  loops  is  vectorized.  In  the 
first,  the  interaction  contributions  of  panicle  /  are  summed  to  panicle  i.  The  second  updates 
particle  j  by  contributions  from  particle  i. 

The  particle  interactions  themselves  are  ordered  through  a  linked-list  [17]  which  efficiently 
determines  the  neighbor  particles  contributing  to  the  forces  on  each  particle.  Only  those 
neighbor  particles  within  the  compact  support  of  the  smoothing  function  need  to  be 

considered.  The  n2  interactions  that  result  from  direct  application  of  the  SPH  formalism, 
without  consideration  of  the  finite  support  of  the  smoothing  function,  are  reduced  to  order 
n  log  n  interactions  by  means  of  the  linked-list.  The  interaction  lists  are  further  reduced  in 
length  by  taking  advantage  of  the  symmetry  in  the  interactions  and  the  activity  flags. 
Reflective  boundaries  are  incorporated  through  the  linked-list  routines  by  means  of  “ghost” 
particles,  which  are  fictitious  particles  introduced  just  outside  the  computational  domain  to 
balance  the  forces  on  boundary  particles  and  mimic  the  effect  of  perfectly  reflecting 
boundaries.  Outflow  boundaries  have  also  been  implemented. 

The  basic  linked  list  algorithm  is  composed  of  two  separate  routines.  The  first  routine 
performs  several  book-keeping  functions  and  is  executed  only  once  per  problem  cycle  (time 
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step).  This  routine  begins  by  calculating  a  box  number  for  each  particle  based  on  its 
position  and  the  position  of  a  regular  grid  overlaid  on  top  of  the  computational  domain  (the 
area  or  volume  that  bound  the  computational  perimeter)  if  the  particle  is  within  twice  the 
smoothing  length  (2/t)  of  the  computational  boundary.  Finally,  a  box-ordered  linked-list  is 
assembled  containing  all  particles  (real  and  ghost)  in  order  of  increasing  box  number.  The 
number  of  particles  in  each  box  and  the  box  offset  (first  position)  in  the  linked-list  are  also 
stored.  The  box-ordered  linked-list,  box  offsets,  and  number  of  particles  per  box  are  used 
in  the  second  routine  described  below. 

The  second  subroutine  is  used  to  find  all  nearest  neighbors  for  particle  i  in  the  hydrodynamic 
calculation  loop  which  are  subsequently  returned  in  a  nearest  neighbor  linked-list.  This  is 
accomplished  by  looping  over  all  particles  contained  in  the  adjacent  boxes  (defined  in  the 
box-ordered  linked-list)  that  surround  the  box  containing  particle  i.  If  the  nearest  neighbor 
index,;',  for  one  of  these  particles  is  greater  than  i  symmetry),  and  either  particle  i  or  j 
is  active,  then  particle  j  is  tested  for  interaction  proximity  to  particle  i.  If  all  of  these 

conditions  are  true  ( j>i,  i  or  j  active,  and  |  r,  -  r,  |  <  2/t),  then  particle  j  is  added  to  the 
nearest  neighbor  linked-list  of  particle  i. 

CALCULATIONS 

The  Noh  Problem  -  The  uniform  implosion  of  an  ideal  (7  =  5/3)  gas  was  conceived  by  Noh  [  10] 
as  a  stringent  test  problem  for  shock  codes.  Initially,  the  gas  is  moving  radially  inward  at  unit 
speed,  unit  density  and  zero  internal  energy.  Noh  found  the  analytic  solution  to  be  a  shock 
moving  radially  outward  at  speed  1/3.  In  spherical  geometry  the  gas  behind  the  shock  has 
particle  speed  0,  specific  internal  energy  1/2  and  density  64.  The  value  of  64  is  due  to  a 
16-fold  increase  from  adiabatic  compression  and  a  4-fold  increase  across  a  strong  shock  for  a 
monatomic  gas.  Our  calculation  used  one-eighth  of  a  sphere  in  three-dimensional  Cartesian 
coordinates  and  three  reflecting  planes.  Particles  were  placed  within  this  domain  in  a  regula. 
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cubic  array  and  then  randomly  perturbed  with  maximum  excursion  oih/H.  The  initial  radius 
of  the  particle  cloud  was  70.  The  smoothing  length  was  set  to  1  with  1  particle-per-h  in  each 
coordinate  direction  giving  roughly  200,000  particles,  including  ghost  particles.  Each  particle 
was  given  unit  density,  unit  speed  inward  and  zero  internal  energy.  Results  of  the  SPH 
calculation  are  shown  in  Figure  la  where  the  density  is  plotted  as  a  function  of  radius  for  each 
particle  at  time  =  48.  Notice  that  all  SPH  particles  fall  on  one  curve  showing  that  perfect 
symmetry  is  achieved  in  the  calculation.  This  is  the  result  to  be  expected  as  there  is  no  spatial 
mesh  which  can  bias  the  solution  along  gridlines.  The  shock  is  in  the  right  place  and  the 
density  dip  (energy  spike)  near  the  origin  is  kept  small  by  the  Noh  wall-heating  term.  The 
calculation  took  8  hours  to  run  on  a  CRAY2  machine.  This  is  a  relatively  long  time,  we 
suspect,  compared  to  other  methods.  The  reason  for  the  slowness  is  due  to  the  implosion 
nature  of  the  problem  coupled  with  our  linked— list  neighbor  algorithm.  As  the  gas  continues 
to  move  radially  inward,  the  calculational  time  increases  dramatically  as  the  number  of 
interacting  neighbors  for  each  particle  increases  from  32  to  2000  because  the  particles  are 
piling  up  near  the  origin  but  the  smoothing  length  remains  fixed.  An  SPH  calculation  with 
variable  smoothing  length  [18]  would  prove  much  more  efficient  for  this  problem.  For 
explosions  and  rarefactions  the  variable  smoothing  length  is  often  required  in  order  to 
maintain  resolution  in  an  expanding  particle  cloud.  We  also  present  results  of  the 
“cylindrical”  Noh  problem  in  Figure  lb.  In  this  geometry  the  gas  behind  the  shock  has 
particle  speed  0,  specific  internal  energy  1/2  and  density  16.  The  calculation  was  run  to 
time  =  60  in  two-dimensional  Cartesian  coordinates  and  required  10,000  particles  and  30 
minutes  of  CRAY2  time. 

Cylinder  Impact  Test  -  Numerical  simulation  of  the  deformation  of  a  metal  cylinder  resulting 
from  normal  impact  against  a  flat,  rigid  surface  is  often  used  to  test  constitutive  models  in 
codes.  There  is  ample  experimental  data  and  the  tests  are  simple  yet  stringent.  We  have 
modeled  an  ARMCO  Iron  cylinder  with  speed  221  m/s  impacting  a  perfectly  reflecting 
surface  using  SPH.  One-quarter  of  the  cylinder  and  two  reflecting  planes  were  used  to 
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a.  3D  spherical  implosion  at  time=48. 


b.  2D  cylindrical  implosion  at  time=60. 

Figure  1.  Density  profile  for  the  Noh  implosion  problem  as  a  function  of  radius 
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model  the  cylinder.  A  third  reflecting  plane  represented  a  perfectly  rigid  boundary.  The 
initial  length  of  the  Iron  rod  was  Lo= 2.54  cm  and  the  initial  diameter  was  Z)o=0.76  cm.  The 
smoothing  length  was  chosen  to  be  h= 0.076  cm  with  2  particles-per-h  in  each  coordinate 
direction.  A  total  of  24,455  particles  were  used.  The  yield  strength  (Yo)  and  shear  modulus 
(p)  of  the  Iron  were  taken  to  be  6  and  0.1  Kb  respectively.  An  initial  density  of  7.89  g/cc  was 
used.  Figure  2  shows  the  final  shape  of  the  computed  cylinder  (2a)  next  to  a  photograph  of  the 
experimental  [19j  article  (2b).  The  the  ratio  of  final  length  and  initial  length  of  the  actual 
experimental  rod  was  LelLo  =  0.78.  The  calculation  gave  LcILo  =  0.79.  The  diameter  ratios 
were  De/Do  =  1.80  for  the  experiment  and  Dc/Do  — 1.55  for  the  calculation,  showing  that  the 
simulation  has  underestimated  the  bulge  near  the  base  of  the  rod.  The  calculation  required  3 
hours  of  CRAY2  time.  This  relatively  long  run  time  is  due  to  the  small  impact  speed  (0.221 
km/s)  compared  to  sound  speed  in  iron  (4.0  km/s)  which  controls  the  time  step. 

Hypervelocity  Impact  -  Figure  3a  shows  the  SPH  calculated  debris  cloud  resulting  from  the 
normal  impact  of  a  3  g  Copper  disk  (11.18  mm  dia  x  3.45  mm  thick)  on  a  2.87  mm  thick 
Aluminum  bumper  plate  at  5.55  km/s.  Figure  3b  is  a  radiograph  of  the  actual  cloud  taken 
from  Piekutowski  [20].  The  experimental  impact  was  not  exactly  normal,  the  Copper  disk 
having  a  5.4  deg  yaw.  We  took  the  smoothing  length  to  be  h= 0.20  mm  and  2  particles  per  h 
giving  10,000  particles  total  in  the  calculation.  A  Gruneisen  equation  of  state  with  Copper 
Hugoniot  Us=0.39+1 .5Up  and  Aluminum  Hugoniot  Us=0.53+I ,5Up  was  used  to  describe  the 
lead  in  compression.  The  shear  modulus  (p)  and  yield  strength  (Yo)  for  Copper  was  taken  to 
be  0.46  Mb  and  4.50  Kb  respectively.  For  Aluminum  we  used  p  =  0.25  Mb  and  Yo=5.50  Kb. 
The  calculation  took  900  cycles  and  1 .8  c.p.u.  hours  on  a  CRAY2.  The  peculiar  shape  of  the 
Aluminum  debris  cloud  is  captured  by  the  simulation.  Figures  3c  and  3d  show  three- 
dimensional  views  of  the  particles.  Figure  3c  shows  only  Aluminum  bumper  plate  particles 
and  3d  shows  only  Copper  projectile  particles. 
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Experimental  photograph  [19]. 
Figure  2.  Cylinder  impact  test. 
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c.  Aluminum  particles  in  3D  perspective 
(  calculation) 


d.  Copper  particles  in  3D  perspective 
(calculation) 


Figure  3.  Debris  cloud  for  the  impact  of  a  Copper  pellet  on  an  Aluminum  bumper  plate. 
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DISCUSSION 


The  three-dimensional  Smoothed  Particle  Hydrodynamics  code  MAGI  has  been  described 
and  three  calculations  presented.  Results  of  the  calculations  are  in  reasonably  good 
•  agreement  with  experiment  and  analysis  showing  that  SPH  can  be  applied  to  low  speed 

■  impacts  as  well  as  hypervelocity  collisions  where  material  strength  is  unimportant. 

Advantages  of  the  method  are  its  robustness,  conceptual  simplicity,  ease  of  adding  new 
physics,  a  natural  treatment  of  void  and  the  ability  to  handle  high  strains  in  a  pure 
Lagrangian  frame.  Tracking  of  debris  clouds  resulting  from  hypervelocity  impacts  is  a 
particularly  important  advantage  of  the  method.  The  run  times  appear  to  be  larger  than  for 
Eulerian  codes  although  no  direct  comparisons  have  been  made.  A  variable  smoothing 
length  formulation  of  SPH  would  dramatically  improve  the  running  time  for  the  Noh 
implosion  and  the  lead  impact  problems  presented  here.  More  fundamentally,  SPH  appears 
to  be  readily  parallelizable.  If  so,  a  one  or  two  ouer  of  magnitude  speed  up  is  possible  on 
todays  massively  parallel  machines.  This  is  an  important  area  of  research. 
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